from casetoolbox.casestab import casestab
from casetoolbox.casestab import math_functions as mf
import os
import numpy as np
from matplotlib import pyplot as plt
# Get current work directory
work_folder = os.getcwd()
# Set result folder
res_folder = 'SteadyState/results_6x6'
# Setup rotor models for each operational point
rotor_models = casestab.rotor_models(os.path.join(work_folder,'DTU10MW_6x6.json'))
# Compute the steady for all operational points
rotor_models.steady_state_computation()
# Save the steady for all operational points
rotor_models.save_steady_state_results('/' + res_folder + '/ind/')
# Copy for shorter notation in plotting below
ops = rotor_models.ops
pwr = rotor_models.pwr
thr = rotor_models.thr
models = rotor_models.models
# Plot comparison to HS2 for each operational point
for iops in range(len(models)):
    # Generate output data
    defl_data = models[iops].substructures[1].create_data_for_deflection_state()
    bem_data = models[iops].rotors[0].create_data_for_BEM_results()
    # Wind speed for labels
    wsp = ops[iops,0]
    # Read HS2 results
    hs2_ind = np.loadtxt(os.path.join(work_folder,'HAWCStab2_results/dtu10mw_u{:d}.ind'.format(int(1000*wsp))),skiprows = 1)
    # Plot deflections and forces in the rotor plane
    fig1=plt.figure(figsize=[12,6])
    axs1=fig1.subplots(3,2)
    axs1_2 = []
    for i in range(3):
        axdummy = []
        for j in range(2):
            axdummy.append(axs1[i,j].twinx())
        axs1_2.append(axdummy)
    axs1[0,0].plot(defl_data[:,0],defl_data[:,4],'.-',hs2_ind[:,0],hs2_ind[:,9],'.-')
    axs1[1,0].plot(defl_data[:,0],defl_data[:,5],'.-',hs2_ind[:,0],hs2_ind[:,10],'.-')
    axs1[2,0].plot(defl_data[:,0],defl_data[:,6],'.-',hs2_ind[:,0],hs2_ind[:,11],'.-')
    axs1[0,1].plot(defl_data[:,0],defl_data[:,9],'.-',hs2_ind[:,0],np.degrees(hs2_ind[:,34]*hs2_ind[:,37]),'.-')
    axs1[1,1].plot(bem_data[:,0],bem_data[:,5]*1e-3,'.-',hs2_ind[:,0],1e-3*hs2_ind[:,6],'.-')
    axs1[2,1].plot(bem_data[:,0],bem_data[:,6]*1e-3,'.-',hs2_ind[:,0],1e-3*hs2_ind[:,7],'.-')
    axs1[0,0].set_ylabel('x deflection [m]')
    axs1[1,0].set_ylabel('y deflection [m]')
    axs1[2,0].set_ylabel('z deflection [m]')
    axs1[0,1].set_ylabel('z rotation [deg]')
    axs1[1,1].set_ylabel('Inplane force [kN/m]')
    axs1[2,1].set_ylabel('Out-of-plane force [kN/m]')
    axs1[0,0].legend(['CS','HS2'])
    axs1_2[0][0].plot(bem_data[1:-1,0],1e3*(np.interp(bem_data[1:-1,0],defl_data[:,0],defl_data[:,4])-hs2_ind[:,9]),'.-g')
    axs1_2[1][0].plot(bem_data[1:-1,0],1e3*(np.interp(bem_data[1:-1,0],defl_data[:,0],defl_data[:,5])-hs2_ind[:,10]),'.-g')
    axs1_2[2][0].plot(bem_data[1:-1,0],1e3*(np.interp(bem_data[1:-1,0],defl_data[:,0],defl_data[:,6])-hs2_ind[:,11]),'.-g')
    axs1_2[0][1].plot(bem_data[1:-1,0],np.interp(bem_data[1:-1,0],defl_data[:,0],defl_data[:,9])-np.degrees(hs2_ind[:,34]*hs2_ind[:,37]),'.-g')
    axs1_2[1][1].plot(bem_data[1:-1,0],bem_data[1:-1,5]-hs2_ind[:,6],'.-g')
    axs1_2[2][1].plot(bem_data[1:-1,0],bem_data[1:-1,6]-hs2_ind[:,7],'.-g')
    axs1_2[0][0].set_ylabel('Diff CS - HS2 [mm]')
    axs1_2[1][0].set_ylabel('Diff CS - HS2 [mm]')
    axs1_2[2][0].set_ylabel('Diff CS - HS2 [mm]')
    axs1_2[0][1].set_ylabel('Diff CS - HS2 [deg]')
    axs1_2[1][1].set_ylabel('Diff CS - HS2 [N/m]')
    axs1_2[2][1].set_ylabel('Diff CS - HS2 [N/m]')
    axs1[2,0].set_xlabel('Spanwise coordinate [m]')
    axs1[2,1].set_xlabel('Spanwise coordinate [m]')
    plt.tight_layout()
    # plt.show()
    fig1.savefig(os.path.join(work_folder, res_folder + '/deflforce{:d}ms.png'.format(int(wsp))),dpi=300)

    plt.close('all')

    veltri = models[iops].rotors[0].blades[0].states_and_forces()

    # Plot velocity triangles
    fig2=plt.figure(figsize=[12,6])
    axs2=fig2.subplots(3,2)
    axs2_2 = []
    for i in range(3):
        axdummy = []
        for j in range(2):
            axdummy.append(axs2[i,j].twinx())
        axs2_2.append(axdummy)
    axs2[0,0].plot(bem_data[:,0],bem_data[:,8],'.-', hs2_ind[:,0],hs2_ind[:, 16],'.-')
    axs2[1,0].plot(bem_data[:,0],bem_data[:,9],'.-', hs2_ind[:,0],hs2_ind[:, 17],'.-')
    axs2[2,0].plot(bem_data[:,0],bem_data[:,3],'.-' ,hs2_ind[:,0],np.degrees(hs2_ind[:, 4]),'.-')
    axs2[0,1].plot(bem_data[:,0],bem_data[:,4],'.-' ,hs2_ind[:,0],hs2_ind[:, 5],'.-')
    axs2[1,1].plot(bem_data[:,0],bem_data[:,1],'.-' ,hs2_ind[:,0],hs2_ind[:, 1],'.-')
    axs2[2,1].plot(bem_data[:,0],bem_data[:,2],'.-' ,hs2_ind[:,0],hs2_ind[:, 2],'.-')
    axs2[0,0].set_ylabel('CL [-]')
    axs2[1,0].set_ylabel('CD [-]')
    axs2[2,0].set_ylabel('AoA [deg]')
    axs2[0,1].set_ylabel('Relative speed [m/s]')
    axs2[1,1].set_ylabel('a [-]')
    axs2[2,1].set_ylabel('ap [-]')
    axs2[0,0].legend(['CS','HS2'])
    axs2_2[0][0].plot(bem_data[1:-1,0],bem_data[1:-1,8]-hs2_ind[:, 16],'.-g')
    axs2_2[1][0].plot(bem_data[1:-1,0],bem_data[1:-1,9]-hs2_ind[:, 17],'.-g')
    axs2_2[2][0].plot(bem_data[1:-1,0],bem_data[1:-1,3]-np.degrees(hs2_ind[:, 4]),'.-g')
    axs2_2[0][1].plot(bem_data[1:-1,0],bem_data[1:-1,4]-hs2_ind[:, 5],'.-g')
    axs2_2[1][1].plot(bem_data[1:-1,0],bem_data[1:-1,1]-hs2_ind[:, 1],'.-g')
    axs2_2[2][1].plot(bem_data[1:-1,0],bem_data[1:-1,2]-hs2_ind[:, 2],'.-g')
    axs2_2[0][0].set_ylabel('Diff CS - HS2 [-]')
    axs2_2[1][0].set_ylabel('Diff CS - HS2 [-]')
    axs2_2[2][0].set_ylabel('Diff CS - HS2 [deg]')
    axs2_2[0][1].set_ylabel('Diff CS - HS2 [m/s]')
    axs2_2[1][1].set_ylabel('Diff CS - HS2 [-]')
    axs2_2[2][1].set_ylabel('Diff CS - HS2 [-]')
    axs2[2,0].set_xlabel('Spanwise position [m]')
    axs2[2,1].set_xlabel('Spanwise position [m]')
    plt.tight_layout()
    # plt.show()
    fig2.savefig(os.path.join(work_folder,res_folder + '/veltri{:d}ms.png'.format(int(wsp))),dpi=300)

    plt.close('all')



hs2_pwr = np.loadtxt(os.path.join(work_folder,'HAWCStab2_results/hs2.pwr'),skiprows = 1)

# Plot power curve
fig = plt.figure()
ax = fig.gca()
ax2 = ax.twinx()

ax.plot(ops[:,0],pwr*1e-3,'o-',hs2_pwr[:,0],hs2_pwr[:,1],'o-')
ax2.plot(ops[:,0],pwr*1e-3-hs2_pwr[:,1],'o-g')

ax.set_ylabel('Power [kW]')
ax.legend(['CS','HS2'])
ax2.set_ylabel('Diff CS - HS2 [kW]')
ax.set_xlabel('Wind speed [m/s]')

ax.set_xlim([ops[0,0],ops[-1,0]])

plt.tight_layout()
#plt.show()
fig.savefig(os.path.join(work_folder,res_folder + '/power.png'),dpi=300)

# Plot thrust curve
fig = plt.figure()
ax = fig.gca()
ax2 = ax.twinx()

ax.plot(ops[:,0],thr*1e-3,'o-',hs2_pwr[:,0],hs2_pwr[:,2],'o-')
ax2.plot(ops[:,0],thr*1e-3-hs2_pwr[:,2],'o-g')

ax.set_ylabel('Thrust [kN]')
ax.legend(['CS','HS2'])
ax2.set_ylabel('Diff CS - HS2 [kN]')
ax.set_xlabel('Wind speed [m/s]')

ax.set_xlim([ops[0,0],ops[-1,0]])

plt.tight_layout()
#plt.show()
fig.savefig(os.path.join(work_folder,res_folder + '/thrust.png'),dpi=300)

# 在现有代码基础上添加以下功能

# =============================================================================
# 增强载荷输出功能
# =============================================================================

def save_detailed_loads_data(rotor_models, work_folder, res_folder):
    """
    保存详细的载荷数据到文件
    """
    import pandas as pd
    
    # 创建载荷数据文件夹
    loads_folder = os.path.join(work_folder, res_folder, 'loads_data')
    os.makedirs(loads_folder, exist_ok=True)
    
    ops = rotor_models.ops
    models = rotor_models.models
    
    # 为每个工况保存详细载荷数据
    for iops in range(len(models)):
        wsp = ops[iops, 0]
        
        # 提取载荷数据
        defl_data = models[iops].substructures[1].create_data_for_deflection_state()
        bem_data = models[iops].rotors[0].create_data_for_BEM_results()
        
        # 创建载荷数据DataFrame
        loads_df = pd.DataFrame({
            'Spanwise_position_[m]': bem_data[:, 0],
            'Inplane_force_[N/m]': bem_data[:, 5],
            'Out_of_plane_force_[N/m]': bem_data[:, 6],
            'Tangential_force_[N/m]': bem_data[:, 5],  # 面内力
            'Normal_force_[N/m]': bem_data[:, 6],        # 面外力
            'Relative_velocity_[m/s]': bem_data[:, 4],
            'Angle_of_attack_[deg]': np.degrees(bem_data[:, 3]),
            'CL': bem_data[:, 8],
            'CD': bem_data[:, 9],
            'Induction_factor_a': bem_data[:, 1],
            'Induction_factor_ap': bem_data[:, 2]
        })
        
        # 保存为CSV文件
        csv_filename = os.path.join(loads_folder, f'loads_wsp_{wsp:.1f}ms.csv')
        loads_df.to_csv(csv_filename, index=False)
        
        # 保存为Excel文件（包含多个工作表）
        excel_filename = os.path.join(loads_folder, f'loads_wsp_{wsp:.1f}ms.xlsx')
        with pd.ExcelWriter(excel_filename, engine='openpyxl') as writer:
            loads_df.to_excel(writer, sheet_name='Loads', index=False)
            
            # 添加变形数据
            defl_df = pd.DataFrame({
                'Spanwise_position_[m]': defl_data[:, 0],
                'x_deflection_[m]': defl_data[:, 4],
                'y_deflection_[m]': defl_data[:, 5],
                'z_deflection_[m]': defl_data[:, 6],
                'x_rotation_[deg]': np.degrees(defl_data[:, 7]),
                'y_rotation_[deg]': np.degrees(defl_data[:, 8]),
                'z_rotation_[deg]': np.degrees(defl_data[:, 9])
            })
            defl_df.to_excel(writer, sheet_name='Deflections', index=False)
            
            # 添加性能数据
            perf_df = pd.DataFrame({
                'Wind_speed_[m/s]': [wsp],
                'Power_[W]': [rotor_models.pwr[iops]],
                'Thrust_[N]': [rotor_models.thr[iops]],
                'Rotor_speed_[rpm]': [models[iops].rotors[0].omega * 30 / np.pi],
                'Tip_speed_ratio': [models[iops].rotors[0].tip_speed_ratio]
            })
            perf_df.to_excel(writer, sheet_name='Performance', index=False)

def plot_enhanced_loads_analysis(rotor_models, work_folder, res_folder):
    """
    绘制增强的载荷分析图表
    """
    ops = rotor_models.ops
    models = rotor_models.models
    
    for iops in range(len(models)):
        wsp = ops[iops, 0]
        defl_data = models[iops].substructures[1].create_data_for_deflection_state()
        bem_data = models[iops].rotors[0].create_data_for_BEM_results()
        
        # 读取HAWCStab2对比数据
        hs2_ind = np.loadtxt(os.path.join(work_folder, 'HAWCStab2_results/dtu10mw_u{:d}.ind'.format(int(1000*wsp))), skiprows=1)
        
        # 1. 载荷分布详细分析
        fig1 = plt.figure(figsize=[16, 12])
        axs1 = fig1.subplots(4, 3)
        
        # 载荷分布图
        axs1[0,0].plot(bem_data[:,0], bem_data[:,5]*1e-3, '.-', label='CASEStab')
        axs1[0,0].plot(hs2_ind[:,0], 1e-3*hs2_ind[:,6], '.-', label='HAWCStab2')
        axs1[0,0].set_ylabel('Inplane Force [kN/m]')
        axs1[0,0].set_title('Inplane Force Distribution')
        axs1[0,0].legend()
        axs1[0,0].grid(True)
        
        axs1[0,1].plot(bem_data[:,0], bem_data[:,6]*1e-3, '.-', label='CASEStab')
        axs1[0,1].plot(hs2_ind[:,0], 1e-3*hs2_ind[:,7], '.-', label='HAWCStab2')
        axs1[0,1].set_ylabel('Out-of-plane Force [kN/m]')
        axs1[0,1].set_title('Out-of-plane Force Distribution')
        axs1[0,1].legend()
        axs1[0,1].grid(True)
        
        # 载荷差异分析
        axs1[0,2].plot(bem_data[1:-1,0], (bem_data[1:-1,5] - hs2_ind[:,6])*1e-3, '.-g')
        axs1[0,2].set_ylabel('Force Difference [kN/m]')
        axs1[0,2].set_title('Inplane Force Difference')
        axs1[0,2].grid(True)
        
        # 2. 气动系数分析
        axs1[1,0].plot(bem_data[:,0], bem_data[:,8], '.-', label='CASEStab')
        axs1[1,0].plot(hs2_ind[:,0], hs2_ind[:,16], '.-', label='HAWCStab2')
        axs1[1,0].set_ylabel('Lift Coefficient [-]')
        axs1[1,0].set_title('Lift Coefficient Distribution')
        axs1[1,0].legend()
        axs1[1,0].grid(True)
        
        axs1[1,1].plot(bem_data[:,0], bem_data[:,9], '.-', label='CASEStab')
        axs1[1,1].plot(hs2_ind[:,0], hs2_ind[:,17], '.-', label='HAWCStab2')
        axs1[1,1].set_ylabel('Drag Coefficient [-]')
        axs1[1,1].set_title('Drag Coefficient Distribution')
        axs1[1,1].legend()
        axs1[1,1].grid(True)
        
        # 攻角分布
        axs1[1,2].plot(bem_data[:,0], np.degrees(bem_data[:,3]), '.-', label='CASEStab')
        axs1[1,2].plot(hs2_ind[:,0], np.degrees(hs2_ind[:,4]), '.-', label='HAWCStab2')
        axs1[1,2].set_ylabel('Angle of Attack [deg]')
        axs1[1,2].set_title('Angle of Attack Distribution')
        axs1[1,2].legend()
        axs1[1,2].grid(True)
        
        # 3. 变形分析
        axs1[2,0].plot(defl_data[:,0], defl_data[:,4], '.-', label='CASEStab')
        axs1[2,0].plot(hs2_ind[:,0], hs2_ind[:,9], '.-', label='HAWCStab2')
        axs1[2,0].set_ylabel('x Deflection [m]')
        axs1[2,0].set_title('x Deflection Distribution')
        axs1[2,0].legend()
        axs1[2,0].grid(True)
        
        axs1[2,1].plot(defl_data[:,0], defl_data[:,5], '.-', label='CASEStab')
        axs1[2,1].plot(hs2_ind[:,0], hs2_ind[:,10], '.-', label='HAWCStab2')
        axs1[2,1].set_ylabel('y Deflection [m]')
        axs1[2,1].set_title('y Deflection Distribution')
        axs1[2,1].legend()
        axs1[2,1].grid(True)
        
        axs1[2,2].plot(defl_data[:,0], defl_data[:,6], '.-', label='CASEStab')
        axs1[2,2].plot(hs2_ind[:,0], hs2_ind[:,11], '.-', label='HAWCStab2')
        axs1[2,2].set_ylabel('z Deflection [m]')
        axs2[2,2].set_title('z Deflection Distribution')
        axs1[2,2].legend()
        axs1[2,2].grid(True)
        
        # 4. 扭转分析
        axs1[3,0].plot(defl_data[:,0], np.degrees(defl_data[:,9]), '.-', label='CASEStab')
        axs1[3,0].plot(hs2_ind[:,0], np.degrees(hs2_ind[:,34]*hs2_ind[:,37]), '.-', label='HAWCStab2')
        axs1[3,0].set_ylabel('z Rotation [deg]')
        axs1[3,0].set_xlabel('Spanwise Position [m]')
        axs1[3,0].set_title('z Rotation Distribution')
        axs1[3,0].legend()
        axs1[3,0].grid(True)
        
        # 相对速度
        axs1[3,1].plot(bem_data[:,0], bem_data[:,4], '.-', label='CASEStab')
        axs1[3,1].plot(hs2_ind[:,0], hs2_ind[:,5], '.-', label='HAWCStab2')
        axs1[3,1].set_ylabel('Relative Velocity [m/s]')
        axs1[3,1].set_xlabel('Spanwise Position [m]')
        axs1[3,1].set_title('Relative Velocity Distribution')
        axs1[3,1].legend()
        axs1[3,1].grid(True)
        
        # 诱导因子
        axs1[3,2].plot(bem_data[:,0], bem_data[:,1], '.-', label='CASEStab')
        axs1[3,2].plot(hs2_ind[:,0], hs2_ind[:,1], '.-', label='HAWCStab2')
        axs1[3,2].set_ylabel('Induction Factor a [-]')
        axs1[3,2].set_xlabel('Spanwise Position [m]')
        axs1[3,2].set_title('Axial Induction Factor')
        axs1[3,2].legend()
        axs1[3,2].grid(True)
        
        plt.tight_layout()
        fig1.savefig(os.path.join(work_folder, res_folder, f'enhanced_loads_analysis_{wsp:.1f}ms.png'), dpi=300)
        plt.close(fig1)

def generate_loads_summary_report(rotor_models, work_folder, res_folder):
    """
    生成载荷汇总报告
    """
    ops = rotor_models.ops
    models = rotor_models.models
    
    # 创建报告文件夹
    report_folder = os.path.join(work_folder, res_folder, 'loads_reports')
    os.makedirs(report_folder, exist_ok=True)
    
    # 汇总数据
    summary_data = []
    
    for iops in range(len(models)):
        wsp = ops[iops, 0]
        bem_data = models[iops].rotors[0].create_data_for_BEM_results()
        
        # 计算载荷统计
        max_inplane_force = np.max(bem_data[:, 5])
        max_outplane_force = np.max(bem_data[:, 6])
        total_thrust = rotor_models.thr[iops]
        total_power = rotor_models.pwr[iops]
        
        # 计算载荷积分（总力）
        span_positions = bem_data[:, 0]
        inplane_forces = bem_data[:, 5]
        outplane_forces = bem_data[:, 6]
        
        # 使用梯形积分计算总力
        total_inplane_force = np.trapz(inplane_forces, span_positions)
        total_outplane_force = np.trapz(outplane_forces, span_positions)
        
        summary_data.append({
            'Wind_Speed_[m/s]': wsp,
            'Max_Inplane_Force_[N/m]': max_inplane_force,
            'Max_Outplane_Force_[N/m]': max_outplane_force,
            'Total_Inplane_Force_[N]': total_inplane_force,
            'Total_Outplane_Force_[N]': total_outplane_force,
            'Total_Thrust_[N]': total_thrust,
            'Total_Power_[W]': total_power,
            'Power_Coefficient': total_power / (0.5 * 1.225 * np.pi * 89.17**2 * wsp**3),
            'Thrust_Coefficient': total_thrust / (0.5 * 1.225 * np.pi * 89.17**2 * wsp**2)
        })
    
    # 保存汇总报告
    import pandas as pd
    summary_df = pd.DataFrame(summary_data)
    summary_df.to_csv(os.path.join(report_folder, 'loads_summary.csv'), index=False)
    summary_df.to_excel(os.path.join(report_folder, 'loads_summary.xlsx'), index=False)
    
    # 生成载荷报告文本
    with open(os.path.join(report_folder, 'loads_report.txt'), 'w') as f:
        f.write("=== 风力涡轮机载荷分析报告 ===\n\n")
        f.write(f"分析时间: {pd.Timestamp.now()}\n")
        f.write(f"分析工况数: {len(models)}\n\n")
        
        f.write("=== 载荷汇总 ===\n")
        f.write(f"最大面内力: {np.max([d['Max_Inplane_Force_[N/m]'] for d in summary_data]):.2f} N/m\n")
        f.write(f"最大面外力: {np.max([d['Max_Outplane_Force_[N/m]'] for d in summary_data]):.2f} N/m\n")
        f.write(f"最大推力: {np.max([d['Total_Thrust_[N]'] for d in summary_data]):.2f} N\n")
        f.write(f"最大功率: {np.max([d['Total_Power_[W]'] for d in summary_data]):.2f} W\n\n")
        
        f.write("=== 各工况详细数据 ===\n")
        for data in summary_data:
            f.write(f"风速 {data['Wind_Speed_[m/s]']:.1f} m/s:\n")
            f.write(f"  最大面内力: {data['Max_Inplane_Force_[N/m]']:.2f} N/m\n")
            f.write(f"  最大面外力: {data['Max_Outplane_Force_[N/m]']:.2f} N/m\n")
            f.write(f"  总推力: {data['Total_Thrust_[N]']:.2f} N\n")
            f.write(f"  总功率: {data['Total_Power_[W]']:.2f} W\n")
            f.write(f"  功率系数: {data['Power_Coefficient']:.3f}\n")
            f.write(f"  推力系数: {data['Thrust_Coefficient']:.3f}\n\n")

def plot_loads_envelope(rotor_models, work_folder, res_folder):
    """
    绘制载荷包络线
    """
    ops = rotor_models.ops
    models = rotor_models.models
    
    # 收集所有工况的载荷数据
    all_inplane_forces = []
    all_outplane_forces = []
    span_positions = None
    
    for iops in range(len(models)):
        bem_data = models[iops].rotors[0].create_data_for_BEM_results()
        if span_positions is None:
            span_positions = bem_data[:, 0]
        all_inplane_forces.append(bem_data[:, 5])
        all_outplane_forces.append(bem_data[:, 6])
    
    all_inplane_forces = np.array(all_inplane_forces)
    all_outplane_forces = np.array(all_outplane_forces)
    
    # 计算包络线
    max_inplane = np.max(all_inplane_forces, axis=0)
    min_inplane = np.min(all_inplane_forces, axis=0)
    max_outplane = np.max(all_outplane_forces, axis=0)
    min_outplane = np.min(all_outplane_forces, axis=0)
    
    # 绘制包络线图
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
    
    # 面内力包络线
    ax1.fill_between(span_positions, min_inplane*1e-3, max_inplane*1e-3, alpha=0.3, label='Load Envelope')
    ax1.plot(span_positions, max_inplane*1e-3, 'r-', linewidth=2, label='Maximum')
    ax1.plot(span_positions, min_inplane*1e-3, 'b-', linewidth=2, label='Minimum')
    ax1.set_ylabel('Inplane Force [kN/m]')
    ax1.set_title('Inplane Force Envelope')
    ax1.legend()
    ax1.grid(True)
    
    # 面外力包络线
    ax2.fill_between(span_positions, min_outplane*1e-3, max_outplane*1e-3, alpha=0.3, label='Load Envelope')
    ax2.plot(span_positions, max_outplane*1e-3, 'r-', linewidth=2, label='Maximum')
    ax2.plot(span_positions, min_outplane*1e-3, 'b-', linewidth=2, label='Minimum')
    ax2.set_ylabel('Out-of-plane Force [kN/m]')
    ax2.set_xlabel('Spanwise Position [m]')
    ax2.set_title('Out-of-plane Force Envelope')
    ax2.legend()
    ax2.grid(True)
    
    plt.tight_layout()
    fig.savefig(os.path.join(work_folder, res_folder, 'loads_envelope.png'), dpi=300)
    plt.close(fig)

# =============================================================================
# 在主程序中调用增强功能
# =============================================================================

# 在现有的稳态计算后添加以下调用：

# 保存详细载荷数据
save_detailed_loads_data(rotor_models, work_folder, res_folder)

# 绘制增强载荷分析
plot_enhanced_loads_analysis(rotor_models, work_folder, res_folder)

# 生成载荷汇总报告
generate_loads_summary_report(rotor_models, work_folder, res_folder)

# 绘制载荷包络线
plot_loads_envelope(rotor_models, work_folder, res_folder)

print("增强载荷输出功能已完成！")
print(f"详细载荷数据已保存到: {os.path.join(work_folder, res_folder, 'loads_data')}")
print(f"载荷报告已保存到: {os.path.join(work_folder, res_folder, 'loads_reports')}")



"""
A. 初始化和计算
# 设置工作目录和结果文件夹
work_folder = os.getcwd()
res_folder = 'SteadyState/results_6x6'

# 加载转子模型
rotor_models = casestab.rotor_models('DTU10MW_6x6.json')

# 计算稳态
rotor_models.steady_state_computation()

# 保存结果
rotor_models.save_steady_state_results('/' + res_folder + '/ind/')

B. 读取结果
# 读取结果
rotor_models = casestab.rotor_models(os.path.join(work_folder,'DTU10MW_6x6.json'))
rotor_models.load_steady_state_results('/' + res_folder + '/ind/')

# 获取操作点数据
ops = rotor_models.ops
pwr = rotor_models.pwr
thr = rotor_models.thr
models = rotor_models.models
3. 可视化分析
3. 可视化分析
A. 变形和力对比图（图1）
x, y, z方向变形：对比CASEStab和HAWCStab2的叶片变形
扭转角度：叶片扭转角对比
面内和面外力：气动载荷分布对比
误差分析：显示两者之间的差异
B. 速度三角形图（图2）
升力系数CL：气动性能参数
阻力系数CD：气动阻力
攻角AoA：叶片攻角分布
相对速度：来流速度分析
诱导因子a, ap：BEM理论参数
C. 性能曲线
功率曲线：功率-风速关系
推力曲线：推力-风速关系
4. 技术特点
A. 双轴绘图
# 创建双y轴图表
axs1_2 = []
for i in range(3):
    axdummy = []
    for j in range(2):
        axdummy.append(axs1[i,j].twinx())
    axs1_2.append(axdummy)
 主轴：显示实际数值
次轴：显示误差差异
B. 数据插值
   # 使用插值确保数据点对齐
np.interp(bem_data[1:-1,0], defl_data[:,0], defl_data[:,4])
C. 多工况分析
for iops in range(len(models)):
    # 对每个工况进行分析
    wfor iops in range(len(models)):
    # 对每个工况进行分析
    wsp = ops[iops,0]  # 风速
    sp = ops[iops,0]  # 风速
5. 输出结果
A. 图表文件
deflforce{:d}ms.png：变形和力对比图
veltri{:d}ms.png：速度三角形图
power.png：功率曲线
thrust.png：推力曲线
B. 数据文件
稳态计算结果保存到 SteadyState/results_6x6/ind/
6. 验证方法
A. 定量对比
直接数值比较
误差计算和可视化
统计差异分析
B. 可视化验证
曲线形状对比
趋势一致性检查
异常值识别
7. 代码结构优势
模块化设计：清晰的功能分离
自动化处理：批量处理多个工况
详细验证：多维度对比分析
高质量输出：300 DPI图表
错误处理：数据插值确保对齐
8. 应用场景
风力涡轮机设计验证
气动弹性分析
性能预测
软件对比验证
工程报告生成
这个脚本是CASEStab工具箱中用于风力涡轮机分析验证的重要工具，通过详细的对比分析确保计算结果的准确性和可靠性。

载荷计算分析
1. 气动载荷计算
# 从BEM数据中提取载荷信息
bem_data = models[iops].rotors[0].create_data_for_BEM_results()

# 面内力 (Inplane force) - 单位：kN/m
axs1[1,1].plot(bem_data[:,0], bem_data[:,5]*1e-3, '.-', 
               hs2_ind[:,0], 1e-3*hs2_ind[:,6], '.-')

# 面外力 (Out-of-plane force) - 单位：kN/m  
axs1[2,1].plot(bem_data[:,0], bem_data[:,6]*1e-3, '.-',
               hs2_ind[:,0], 1e-3*hs2_ind[:,7], '.-')
2. 载荷类型
A. 面内力 (Inplane force)
B. 面外力 (Out-of-plane force)
C. 扭矩 (Torque)
D. 弯矩 (Bending moment)
E. 剪力 (Shear force)
F. 轴向力 (Axial force)
G. 径向力 (Radial force)
H. 切向力 (Tangential force)
I. 法向力 (Normal force)
B. 整体载荷
# 功率载荷
pwr = rotor_models.pwr  # 功率 [W]

# 推力载荷  
thr = rotor_models.thr  # 推力 [N]
3. 载荷计算过程
A. 稳态计算
# 执行稳态计算，包含载荷计算
rotor_models.steady_state_computation()

￼
￼Add to ChatCtrl+L
Ctrl+K to generate a command
Fix terminal error issue
Generate a name for the chat
Create a new wind turbine example
解释 corotbeam.py 的功能
New Chat
实现Krenk和Couturier平衡单元：更详细的解释

让我详细解释Krenk和Couturier平衡单元的实现。首先让我查看相关的预编译函数文件，以获得更完整的理解：
Read
corotbeam_precompiled_functions.py L1-100
现在让我查看equilibrium_beam_element类的具体实现：
Read
corotbeam.py L707-856
现在让我查看一些相关的数学函数，特别是关于Rodrigues旋转的部分：
Read
math_functions.py L1-100
现在让我查看一些具体的实现细节，特别是关于形状函数计算的部分：
Searched
How are the shape functions computed in the equilibrium beam element? (/home/mich/apps/CASEToolBox/casetoolbox/casestab)
基于对代码的深入分析，我现在可以详细解释Krenk和Couturier平衡单元的实现：
Krenk和Couturier平衡单元详细实现解释
1. 理论基础
Krenk和Couturier平衡单元是一种基于柔度矩阵而非刚度矩阵的有限元方法，特别适用于：
变截面梁结构
复杂材料属性分布
大变形分析
气动弹性耦合问题
2. 核心矩阵构建
H矩阵计算（柔度矩阵积分）
￼
# 计算柔度矩阵的积分for n in range(0, self.norder+1, 2):  # 偶数阶    self.H = self.H + l/(n+1.0) * self.C[n]    tmp = self.C[n] @ T1    self.H = self.H + l/(n+3.0) * T1.T @ tmpfor n in range(1, self.norder+1, 2):  # 奇数阶    tmp = self.C[n] @ T1    self.H = self.H + l/(n+2.0) * (tmp.T + tmp)
物理意义：
H矩阵是柔度矩阵沿单元长度的积分
偶数阶项：直接积分柔度矩阵
奇数阶项：考虑截面变化的影响
T1矩阵（几何变换矩阵）
￼
T1 = np.zeros((6,6))T1[4,0] = -l/2.0  # 弯曲-轴向耦合T1[3,1] =  l/2.0  # 弯曲-轴向耦合
作用：描述轴向变形与弯曲变形的耦合关系
P矩阵（投影矩阵）
￼
P = np.zeros((6,6))P[0,4] =  1.0  # 轴向-弯曲投影P[1,3] = -1.0  # 轴向-弯曲投影
作用：将弯曲变形投影到轴向变形
3. 单元刚度矩阵计算
￼
# 计算H矩阵的逆self.Hinv = np.linalg.inv(self.H)# 构建G^T矩阵self.GT = np.concatenate((-np.eye(6)+T1, np.eye(6)+T1)).T# 单元刚度矩阵self.K = self.GT.T @ self.Hinv @ self.GT
数学原理：
从柔度矩阵C通过积分得到H
通过几何变换GT将局部坐标系转换到全局坐标系
最终得到单元刚度矩阵K
4. 形状函数计算
应变函数Ns
￼
# 0阶项self.Ns.append(self.C[0] @ self.Hinv @ self.GT)# 1到norder阶项for n in range(1, self.norder+1):    self.Ns.append(self.C[n] @ T1 @ self.Hinv @ self.GT -                    self.C[n] @ self.Hinv @ self.GT)# norder+1阶项self.Ns.append(self.C[self.norder] @ T1 @ self.Hinv @ self.GT)
形状函数积分
￼
# 第一次积分Ns1[0] = sum(l/2.0 * (-1.0)**n / (n+1) * self.Ns[n] for n in range(norder+2))# 第二次积分Ns2[0] = sum(l/2.0 * (-1.0)**n / (n+1) * Ns1[n] for n in range(norder+3))
物理意义：
第一次积分：从应变得到位移
第二次积分：从位移得到形状函数
5. 共旋转框架实现
单元坐标系计算
￼
def compute_element_triad_and_position(r, T, Q, dT, dQ, ddT, ddQ):    # 单元轴线    rvec = r[3:] - r[0:3]    l = sqrt(innerproduct(rvec, rvec))        # e3: 单元轴线方向    E[:,2] = rvec / l        # e2: 节点法向量的平均    v = T[:,1] + Q[:,1]    a = innerproduct(v, E[:,2])    d = sqrt(2.0 + 2.0*innerproduct(T[:,1], Q[:,1]) - (a*a))    e2tilde = v - a * E[:,2]    E[:,1] = e2tilde / d        # e1: 叉积得到    E[:,0] = crossproduct(E[:,1], E[:,2])
关键特点：
单元坐标系随变形而旋转
分离刚体运动与局部变形
避免大旋转带来的数值问题
6. 局部变形计算
局部节点旋转和伸长
￼
def update_local_nodal_rotations_elongation(l0, l, rvec0, rvec, dvec, E, T, Q):    # 计算局部变形参数    ql = np.zeros(7)        # 轴向伸长    ql[0] = (l - l0) / l0        # 局部旋转（小角度假设）    # 通过Rodrigues参数计算局部旋转    # ...
7. 形状函数评估
多项式形状函数
￼
def shape_function_matrix(self, zeta):    N = np.zeros((6,12))    for n in range(self.norder+4):        N = N + self.N[n] * zeta**n    return N
特点：
支持高阶单元（通过norder参数）
多项式形式便于积分和微分
分离轴向、弯曲和扭转变形
8. 气动弹性耦合
力分布积分
￼
def compute_element_total_force_matrix(self, elem_model):    Tf = cpf.compute_element_total_force_matrix(        self.ninterval, elem_model.l, self.a, self.b, self.w)    return Tf
功能：
将分布力转换为节点力
支持变截面和复杂力分布
与气动分析模块耦合
9. 数值优势
稳定性：共旋转框架避免大旋转问题
精度：基于柔度矩阵，更准确描述材料行为
效率：预编译函数加速计算
灵活性：支持变截面和复杂材料属性
10. 应用场景
风力涡轮机叶片：大变形、气动弹性耦合
直升机旋翼：复杂载荷、高转速
柔性结构：长细比结构、非线性动力学
这种平衡单元方法为复杂柔性结构的精确分析提供了强大的数值工具，特别适用于风力涡轮机等工程应用。
把上面的回答存成一个文件？

我理解您想要将上面的详细解释保存为文件。由于您当前处于ask模式，我无法直接创建文件，但我可以为您提供完整的文件内容，您可以复制并保存。
以下是完整的文档内容，建议保存为 Krenk_Couturier_Equilibrium_Element_Explanation.md：
￼
# Krenk和Couturier平衡单元详细实现解释## 1. 理论基础Krenk和Couturier平衡单元是一种基于**柔度矩阵**而非刚度矩阵的有限元方法，特别适用于：- 变截面梁结构- 复杂材料属性分布- 大变形分析- 气动弹性耦合问题## 2. 核心矩阵构建### H矩阵计算（柔度矩阵积分）```python# 计算柔度矩阵的积分for n in range(0, self.norder+1, 2):  # 偶数阶    self.H = self.H + l/(n+1.0) * self.C[n]    tmp = self.C[n] @ T1    self.H = self.H + l/(n+3.0) * T1.T @ tmpfor n in range(1, self.norder+1, 2):  # 奇数阶    tmp = self.C[n] @ T1    self.H = self.H + l/(n+2.0) * (tmp.T + tmp)```**物理意义**：- `H`矩阵是柔度矩阵沿单元长度的积分- 偶数阶项：直接积分柔度矩阵- 奇数阶项：考虑截面变化的影响### T1矩阵（几何变换矩阵）```pythonT1 = np.zeros((6,6))T1[4,0] = -l/2.0  # 弯曲-轴向耦合T1[3,1] =  l/2.0  # 弯曲-轴向耦合```**作用**：描述轴向变形与弯曲变形的耦合关系### P矩阵（投影矩阵）```pythonP = np.zeros((6,6))P[0,4] =  1.0  # 轴向-弯曲投影P[1,3] = -1.0  # 轴向-弯曲投影```**作用**：将弯曲变形投影到轴向变形## 3. 单元刚度矩阵计算```python# 计算H矩阵的逆self.Hinv = np.linalg.inv(self.H)# 构建G^T矩阵self.GT = np.concatenate((-np.eye(6)+T1, np.eye(6)+T1)).T# 单元刚度矩阵self.K = self.GT.T @ self.Hinv @ self.GT```**数学原理**：- 从柔度矩阵`C`通过积分得到`H`- 通过几何变换`GT`将局部坐标系转换到全局坐标系- 最终得到单元刚度矩阵`K`## 4. 形状函数计算### 应变函数Ns```python# 0阶项self.Ns.append(self.C[0] @ self.Hinv @ self.GT)# 1到norder阶项for n in range(1, self.norder+1):    self.Ns.append(self.C[n] @ T1 @ self.Hinv @ self.GT -                    self.C[n] @ self.Hinv @ self.GT)# norder+1阶项self.Ns.append(self.C[self.norder] @ T1 @ self.Hinv @ self.GT)```### 形状函数积分```python# 第一次积分Ns1[0] = sum(l/2.0 * (-1.0)**n / (n+1) * self.Ns[n] for n in range(norder+2))# 第二次积分Ns2[0] = sum(l/2.0 * (-1.0)**n / (n+1) * Ns1[n] for n in range(norder+3))```**物理意义**：- 第一次积分：从应变得到位移- 第二次积分：从位移得到形状函数## 5. 共旋转框架实现### 单元坐标系计算```pythondef compute_element_triad_and_position(r, T, Q, dT, dQ, ddT, ddQ):    # 单元轴线    rvec = r[3:] - r[0:3]    l = sqrt(innerproduct(rvec, rvec))        # e3: 单元轴线方向    E[:,2] = rvec / l        # e2: 节点法向量的平均    v = T[:,1] + Q[:,1]    a = innerproduct(v, E[:,2])    d = sqrt(2.0 + 2.0*innerproduct(T[:,1], Q[:,1]) - (a*a))    e2tilde = v - a * E[:,2]    E[:,1] = e2tilde / d        # e1: 叉积得到    E[:,0] = crossproduct(E[:,1], E[:,2])```**关键特点**：- 单元坐标系随变形而旋转- 分离刚体运动与局部变形- 避免大旋转带来的数值问题## 6. 局部变形计算### 局部节点旋转和伸长```pythondef update_local_nodal_rotations_elongation(l0, l, rvec0, rvec, dvec, E, T, Q):    # 计算局部变形参数    ql = np.zeros(7)        # 轴向伸长    ql[0] = (l - l0) / l0        # 局部旋转（小角度假设）    # 通过Rodrigues参数计算局部旋转    # ...```## 7. 形状函数评估### 多项式形状函数```pythondef shape_function_matrix(self, zeta):    N = np.zeros((6,12))    for n in range(self.norder+4):        N = N + self.N[n] * zeta**n    return N```**特点**：- 支持高阶单元（通过`norder`参数）- 多项式形式便于积分和微分- 分离轴向、弯曲和扭转变形## 8. 气动弹性耦合### 力分布积分```pythondef compute_element_total_force_matrix(self, elem_model):    Tf = cpf.compute_element_total_force_matrix(        self.ninterval, elem_model.l, self.a, self.b, self.w)    return Tf```**功能**：- 将分布力转换为节点力- 支持变截面和复杂力分布- 与气动分析模块耦合## 9. 数值优势1. **稳定性**：共旋转框架避免大旋转问题2. **精度**：基于柔度矩阵，更准确描述材料行为3. **效率**：预编译函数加速计算4. **灵活性**：支持变截面和复杂材料属性## 10. 应用场景- **风力涡轮机叶片**：大变形、气动弹性耦合- **直升机旋翼**：复杂载荷、高转速- **柔性结构**：长细比结构、非线性动力学## 11. 关键类和方法### equilibrium_beam_element类- `__init__(l, C, iner_pars)`: 初始化单元- `shape_function_matrix(zeta)`: 计算形状函数- `local_shape_function_matrix(zeta)`: 计算局部形状函数### corotbeam_element_kinematics类- `update_nodal_triads_and_position(q)`: 更新节点位置和方向- `update_local_nodal_rotations_elongation()`: 计算局部变形- `compute_element_triad_and_position()`: 计算单元坐标系### corotbeam_substructure类- `update_elastic_internal_forces_and_stiffness()`: 计算内力和刚度- `update_inertia()`: 更新惯性- `initiate_aeroelastic_coupling()`: 初始化气动弹性耦合## 12. 数学公式总结### 单元刚度矩阵
物理意义：
H矩阵是柔度矩阵沿单元长度的积分
偶数阶项：直接积分柔度矩阵
奇数阶项：考虑截面变化的影响
T1矩阵（几何变换矩阵）
￼
### 形状函数
作用：描述轴向变形与弯曲变形的耦合关系
P矩阵（投影矩阵）
￼
### 局部变形
作用：将弯曲变形投影到轴向变形
3. 单元刚度矩阵计算
￼
这种平衡单元方法为复杂柔性结构的精确分析提供了强大的数值工具，特别适用于风力涡轮机等工程应用。---*文档生成时间：2024年**基于CASEToolBox/CASEStab项目分析*
数学原理：
从柔度矩阵C通过积分得到H
通过几何变换GT将局部坐标系转换到全局坐标系
最终得到单元刚度矩阵K
4. 形状函数计算
应变函数Ns
￼
# 0阶项self.Ns.append(self.C[0] @ self.Hinv @ self.GT)# 1到norder阶项for n in range(1, self.norder+1):    self.Ns.append(self.C[n] @ T1 @ self.Hinv @ self.GT -                    self.C[n] @ self.Hinv @ self.GT)# norder+1阶项self.Ns.append(self.C[self.norder] @ T1 @ self.Hinv @ self.GT)
形状函数积分
￼
# 第一次积分Ns1[0] = sum(l/2.0 * (-1.0)**n / (n+1) * self.Ns[n] for n in range(norder+2))# 第二次积分Ns2[0] = sum(l/2.0 * (-1.0)**n / (n+1) * Ns1[n] for n in range(norder+3))
物理意义：
第一次积分：从应变得到位移
第二次积分：从位移得到形状函数
5. 共旋转框架实现
单元坐标系计算
￼
def compute_element_triad_and_position(r, T, Q, dT, dQ, ddT, ddQ):    # 单元轴线    rvec = r[3:] - r[0:3]    l = sqrt(innerproduct(rvec, rvec))        # e3: 单元轴线方向    E[:,2] = rvec / l        # e2: 节点法向量的平均    v = T[:,1] + Q[:,1]    a = innerproduct(v, E[:,2])    d = sqrt(2.0 + 2.0*innerproduct(T[:,1], Q[:,1]) - (a*a))    e2tilde = v - a * E[:,2]    E[:,1] = e2tilde / d        # e1: 叉积得到    E[:,0] = crossproduct(E[:,1], E[:,2])
关键特点：
单元坐标系随变形而旋转
分离刚体运动与局部变形
避免大旋转带来的数值问题
6. 局部变形计算
局部节点旋转和伸长
￼
def update_local_nodal_rotations_elongation(l0, l, rvec0, rvec, dvec, E, T, Q):    # 计算局部变形参数    ql = np.zeros(7)        # 轴向伸长    ql[0] = (l - l0) / l0        # 局部旋转（小角度假设）    # 通过Rodrigues参数计算局部旋转    # ...
7. 形状函数评估
多项式形状函数
￼
def shape_function_matrix(self, zeta):    N = np.zeros((6,12))    for n in range(self.norder+4):        N = N + self.N[n] * zeta**n    return N
特点：
支持高阶单元（通过norder参数）
多项式形式便于积分和微分
分离轴向、弯曲和扭转变形
8. 气动弹性耦合
力分布积分
￼
def compute_element_total_force_matrix(self, elem_model):    Tf = cpf.compute_element_total_force_matrix(        self.ninterval, elem_model.l, self.a, self.b, self.w)    return Tf
功能：
将分布力转换为节点力
支持变截面和复杂力分布
与气动分析模块耦合
9. 数值优势
稳定性：共旋转框架避免大旋转问题
精度：基于柔度矩阵，更准确描述材料行为
效率：预编译函数加速计算
灵活性：支持变截面和复杂材料属性
10. 应用场景
风力涡轮机叶片：大变形、气动弹性耦合
直升机旋翼：复杂载荷、高转速
柔性结构：长细比结构、非线性动力学
11. 关键类和方法
equilibrium_beam_element类
__init__(l, C, iner_pars): 初始化单元
shape_function_matrix(zeta): 计算形状函数
local_shape_function_matrix(zeta): 计算局部形状函数
corotbeam_element_kinematics类
update_nodal_triads_and_position(q): 更新节点位置和方向
update_local_nodal_rotations_elongation(): 计算局部变形
compute_element_triad_and_position(): 计算单元坐标系
corotbeam_substructure类
update_elastic_internal_forces_and_stiffness(): 计算内力和刚度
update_inertia(): 更新惯性
initiate_aeroelastic_coupling(): 初始化气动弹性耦合
12. 数学公式总结
单元刚度矩阵
￼
K = G^T * H^(-1) * G
形状函数
￼
N(ζ) = Σ(n=0 to norder+3) N[n] * ζ^n
局部变形
￼
u_local = N(ζ) * q_local + 刚体运动
这种平衡单元方法为复杂柔性结构的精确分析提供了强大的数值工具，特别适用于风力涡轮机等工程应用。
文档生成时间：2024年
基于CASEToolBox/CASEStab项目分析
￼
您可以将上述内容复制并保存为 `.md` 文件。如果您需要我帮助创建其他格式的文件（如 `.txt` 或 `.py` 文档字符串），请告诉我。
有很多次的copy，能不能做成一个copy

我理解您的意思，您希望将整个文档内容作为一个完整的代码块，这样您只需要复制一次。以下是完整的文档内容：
￼
# Krenk和Couturier平衡单元详细实现解释## 1. 理论基础Krenk和Couturier平衡单元是一种基于**柔度矩阵**而非刚度矩阵的有限元方法，特别适用于：- 变截面梁结构- 复杂材料属性分布- 大变形分析- 气动弹性耦合问题## 2. 核心矩阵构建### H矩阵计算（柔度矩阵积分）```python# 计算柔度矩阵的积分for n in range(0, self.norder+1, 2):  # 偶数阶    self.H = self.H + l/(n+1.0) * self.C[n]    tmp = self.C[n] @ T1    self.H = self.H + l/(n+3.0) * T1.T @ tmpfor n in range(1, self.norder+1, 2):  # 奇数阶    tmp = self.C[n] @ T1    self.H = self.H + l/(n+2.0) * (tmp.T + tmp)```**物理意义**：- `H`矩阵是柔度矩阵沿单元长度的积分- 偶数阶项：直接积分柔度矩阵- 奇数阶项：考虑截面变化的影响### T1矩阵（几何变换矩阵）```pythonT1 = np.zeros((6,6))T1[4,0] = -l/2.0  # 弯曲-轴向耦合T1[3,1] =  l/2.0  # 弯曲-轴向耦合```**作用**：描述轴向变形与弯曲变形的耦合关系### P矩阵（投影矩阵）```pythonP = np.zeros((6,6))P[0,4] =  1.0  # 轴向-弯曲投影P[1,3] = -1.0  # 轴向-弯曲投影```**作用**：将弯曲变形投影到轴向变形## 3. 单元刚度矩阵计算```python# 计算H矩阵的逆self.Hinv = np.linalg.inv(self.H)# 构建G^T矩阵self.GT = np.concatenate((-np.eye(6)+T1, np.eye(6)+T1)).T# 单元刚度矩阵self.K = self.GT.T @ self.Hinv @ self.GT```**数学原理**：- 从柔度矩阵`C`通过积分得到`H`- 通过几何变换`GT`将局部坐标系转换到全局坐标系- 最终得到单元刚度矩阵`K`## 4. 形状函数计算### 应变函数Ns```python# 0阶项self.Ns.append(self.C[0] @ self.Hinv @ self.GT)# 1到norder阶项for n in range(1, self.norder+1):    self.Ns.append(self.C[n] @ T1 @ self.Hinv @ self.GT -                    self.C[n] @ self.Hinv @ self.GT)# norder+1阶项self.Ns.append(self.C[self.norder] @ T1 @ self.Hinv @ self.GT)```### 形状函数积分```python# 第一次积分Ns1[0] = sum(l/2.0 * (-1.0)**n / (n+1) * self.Ns[n] for n in range(norder+2))# 第二次积分Ns2[0] = sum(l/2.0 * (-1.0)**n / (n+1) * Ns1[n] for n in range(norder+3))```**物理意义**：- 第一次积分：从应变得到位移- 第二次积分：从位移得到形状函数## 5. 共旋转框架实现### 单元坐标系计算```pythondef compute_element_triad_and_position(r, T, Q, dT, dQ, ddT, ddQ):    # 单元轴线    rvec = r[3:] - r[0:3]    l = sqrt(innerproduct(rvec, rvec))        # e3: 单元轴线方向    E[:,2] = rvec / l        # e2: 节点法向量的平均    v = T[:,1] + Q[:,1]    a = innerproduct(v, E[:,2])    d = sqrt(2.0 + 2.0*innerproduct(T[:,1], Q[:,1]) - (a*a))    e2tilde = v - a * E[:,2]    E[:,1] = e2tilde / d        # e1: 叉积得到    E[:,0] = crossproduct(E[:,1], E[:,2])```**关键特点**：- 单元坐标系随变形而旋转- 分离刚体运动与局部变形- 避免大旋转带来的数值问题## 6. 局部变形计算### 局部节点旋转和伸长```pythondef update_local_nodal_rotations_elongation(l0, l, rvec0, rvec, dvec, E, T, Q):    # 计算局部变形参数    ql = np.zeros(7)        # 轴向伸长    ql[0] = (l - l0) / l0        # 局部旋转（小角度假设）    # 通过Rodrigues参数计算局部旋转    # ...```## 7. 形状函数评估### 多项式形状函数```pythondef shape_function_matrix(self, zeta):    N = np.zeros((6,12))    for n in range(self.norder+4):        N = N + self.N[n] * zeta**n    return N```**特点**：- 支持高阶单元（通过`norder`参数）- 多项式形式便于积分和微分- 分离轴向、弯曲和扭转变形## 8. 气动弹性耦合### 力分布积分```pythondef compute_element_total_force_matrix(self, elem_model):    Tf = cpf.compute_element_total_force_matrix(        self.ninterval, elem_model.l, self.a, self.b, self.w)    return Tf```**功能**：- 将分布力转换为节点力- 支持变截面和复杂力分布- 与气动分析模块耦合## 9. 数值优势1. **稳定性**：共旋转框架避免大旋转问题2. **精度**：基于柔度矩阵，更准确描述材料行为3. **效率**：预编译函数加速计算4. **灵活性**：支持变截面和复杂材料属性## 10. 应用场景- **风力涡轮机叶片**：大变形、气动弹性耦合- **直升机旋翼**：复杂载荷、高转速- **柔性结构**：长细比结构、非线性动力学## 11. 关键类和方法### equilibrium_beam_element类- `__init__(l, C, iner_pars)`: 初始化单元- `shape_function_matrix(zeta)`: 计算形状函数- `local_shape_function_matrix(zeta)`: 计算局部形状函数### corotbeam_element_kinematics类- `update_nodal_triads_and_position(q)`: 更新节点位置和方向- `update_local_nodal_rotations_elongation()`: 计算局部变形- `compute_element_triad_and_position()`: 计算单元坐标系### corotbeam_substructure类- `update_elastic_internal_forces_and_stiffness()`: 计算内力和刚度- `update_inertia()`: 更新惯性- `initiate_aeroelastic_coupling()`: 初始化气动弹性耦合## 12. 数学公式总结### 单元刚度矩阵
物理意义：
H矩阵是柔度矩阵沿单元长度的积分
偶数阶项：直接积分柔度矩阵
奇数阶项：考虑截面变化的影响
T1矩阵（几何变换矩阵）
￼
### 形状函数
作用：描述轴向变形与弯曲变形的耦合关系
P矩阵（投影矩阵）
￼
### 局部变形
作用：将弯曲变形投影到轴向变形
3. 单元刚度矩阵计算
￼
这种平衡单元方法为复杂柔性结构的精确分析提供了强大的数值工具，特别适用于风力涡轮机等工程应用。---*文档生成时间：2024年**基于CASEToolBox/CASEStab项目分析*
数学原理：
从柔度矩阵C通过积分得到H
通过几何变换GT将局部坐标系转换到全局坐标系
最终得到单元刚度矩阵K
4. 形状函数计算
应变函数Ns
￼
# 0阶项self.Ns.append(self.C[0] @ self.Hinv @ self.GT)# 1到norder阶项for n in range(1, self.norder+1):    self.Ns.append(self.C[n] @ T1 @ self.Hinv @ self.GT -                    self.C[n] @ self.Hinv @ self.GT)# norder+1阶项self.Ns.append(self.C[self.norder] @ T1 @ self.Hinv @ self.GT)
形状函数积分
￼
# 第一次积分Ns1[0] = sum(l/2.0 * (-1.0)**n / (n+1) * self.Ns[n] for n in range(norder+2))# 第二次积分Ns2[0] = sum(l/2.0 * (-1.0)**n / (n+1) * Ns1[n] for n in range(norder+3))
物理意义：
第一次积分：从应变得到位移
第二次积分：从位移得到形状函数
5. 共旋转框架实现
单元坐标系计算
￼
def compute_element_triad_and_position(r, T, Q, dT, dQ, ddT, ddQ):    # 单元轴线    rvec = r[3:] - r[0:3]    l = sqrt(innerproduct(rvec, rvec))        # e3: 单元轴线方向    E[:,2] = rvec / l        # e2: 节点法向量的平均    v = T[:,1] + Q[:,1]    a = innerproduct(v, E[:,2])    d = sqrt(2.0 + 2.0*innerproduct(T[:,1], Q[:,1]) - (a*a))    e2tilde = v - a * E[:,2]    E[:,1] = e2tilde / d        # e1: 叉积得到    E[:,0] = crossproduct(E[:,1], E[:,2])
关键特点：
单元坐标系随变形而旋转
分离刚体运动与局部变形
避免大旋转带来的数值问题
6. 局部变形计算
局部节点旋转和伸长
￼
def update_local_nodal_rotations_elongation(l0, l, rvec0, rvec, dvec, E, T, Q):    # 计算局部变形参数    ql = np.zeros(7)        # 轴向伸长    ql[0] = (l - l0) / l0        # 局部旋转（小角度假设）    # 通过Rodrigues参数计算局部旋转    # ...
7. 形状函数评估
多项式形状函数
￼
def shape_function_matrix(self, zeta):    N = np.zeros((6,12))    for n in range(self.norder+4):        N = N + self.N[n] * zeta**n    return N
特点：
支持高阶单元（通过norder参数）
多项式形式便于积分和微分
分离轴向、弯曲和扭转变形
8. 气动弹性耦合
力分布积分
￼
def compute_element_total_force_matrix(self, elem_model):    Tf = cpf.compute_element_total_force_matrix(        self.ninterval, elem_model.l, self.a, self.b, self.w)    return Tf
功能：
将分布力转换为节点力
支持变截面和复杂力分布
与气动分析模块耦合
9. 数值优势
稳定性：共旋转框架避免大旋转问题
精度：基于柔度矩阵，更准确描述材料行为
效率：预编译函数加速计算
灵活性：支持变截面和复杂材料属性
10. 应用场景
风力涡轮机叶片：大变形、气动弹性耦合
直升机旋翼：复杂载荷、高转速
柔性结构：长细比结构、非线性动力学
11. 关键类和方法
equilibrium_beam_element类
__init__(l, C, iner_pars): 初始化单元
shape_function_matrix(zeta): 计算形状函数
local_shape_function_matrix(zeta): 计算局部形状函数
corotbeam_element_kinematics类
update_nodal_triads_and_position(q): 更新节点位置和方向
update_local_nodal_rotations_elongation(): 计算局部变形
compute_element_triad_and_position(): 计算单元坐标系
corotbeam_substructure类
update_elastic_internal_forces_and_stiffness(): 计算内力和刚度
update_inertia(): 更新惯性
initiate_aeroelastic_coupling(): 初始化气动弹性耦合
12. 数学公式总结
单元刚度矩阵
￼
K = G^T * H^(-1) * G
形状函数
￼
N(ζ) = Σ(n=0 to norder+3) N[n] * ζ^n
局部变形
￼
u_local = N(ζ) * q_local + 刚体运动
这种平衡单元方法为复杂柔性结构的精确分析提供了强大的数值工具，特别适用于风力涡轮机等工程应用。
文档生成时间：2024年
基于CASEToolBox/CASEStab项目分析
￼
现在您只需要复制上面的整个代码块一次，就可以获得完整的文档内容。
分析一下目前的.py

我来分析一下当前打开的Python文件。让我先读取这个文件的内容：
Read
steady_state_using_HS2_data_translated_6x6.py
现在让我分析这个Python文件的结构和功能：
文件分析：steady_state_using_HS2_data_translated_6x6.py
1. 文件目的
这是一个风力涡轮机稳态分析的脚本，用于：
使用CASEStab计算风力涡轮机的稳态响应
与HAWCStab2 (HS2) 结果进行对比验证
生成详细的对比图表和性能曲线
2. 主要功能模块
A. 初始化和计算
￼
# 设置工作目录和结果文件夹work_folder = os.getcwd()res_folder = 'SteadyState/results_6x6'# 加载转子模型rotor_models = casestab.rotor_models('DTU10MW_6x6.json')# 计算稳态rotor_models.steady_state_computation()# 保存结果rotor_models.save_steady_state_results('/' + res_folder + '/ind/')
B. 数据提取
￼
# 提取变形数据defl_data = models[iops].substructures[1].create_data_for_deflection_state()# 提取BEM（叶片动量理论）数据bem_data = models[iops].rotors[0].create_data_for_BEM_results()# 读取HAWCStab2对比数据hs2_ind = np.loadtxt('HAWCStab2_results/dtu10mw_u{:d}.ind'.format(int(1000*wsp)))
3. 可视化分析
A. 变形和力对比图（图1）
x, y, z方向变形：对比CASEStab和HAWCStab2的叶片变形
扭转角度：叶片扭转角对比
面内和面外力：气动载荷分布对比
误差分析：显示两者之间的差异
B. 速度三角形图（图2）
升力系数CL：气动性能参数
阻力系数CD：气动阻力
攻角AoA：叶片攻角分布
相对速度：来流速度分析
诱导因子a, ap：BEM理论参数
C. 性能曲线
功率曲线：功率-风速关系
推力曲线：推力-风速关系
4. 技术特点
A. 双轴绘图
￼
# 创建双y轴图表axs1_2 = []for i in range(3):    axdummy = []    for j in range(2):        axdummy.append(axs1[i,j].twinx())    axs1_2.append(axdummy)
主轴：显示实际数值
次轴：显示误差差异
B. 数据插值
￼
# 使用插值确保数据点对齐np.interp(bem_data[1:-1,0], defl_data[:,0], defl_data[:,4])
C. 多工况分析
￼
for iops in range(len(models)):    # 对每个工况进行分析    wsp = ops[iops,0]  # 风速
5. 输出结果
A. 图表文件
deflforce{:d}ms.png：变形和力对比图
veltri{:d}ms.png：速度三角形图
power.png：功率曲线
thrust.png：推力曲线
B. 数据文件
稳态计算结果保存到 SteadyState/results_6x6/ind/
6. 验证方法
A. 定量对比
直接数值比较
误差计算和可视化
统计差异分析
B. 可视化验证
曲线形状对比
趋势一致性检查
异常值识别
7. 代码结构优势
模块化设计：清晰的功能分离
自动化处理：批量处理多个工况
详细验证：多维度对比分析
高质量输出：300 DPI图表
错误处理：数据插值确保对齐
8. 应用场景
风力涡轮机设计验证
气动弹性分析
性能预测
软件对比验证
工程报告生成
这个脚本是CASEStab工具箱中用于风力涡轮机分析验证的重要工具，通过详细的对比分析确保计算结果的准确性和可靠性。
有计算载荷吗

让我分析一下这个文件中是否有载荷计算。从代码中我可以看到几个与载荷相关的部分：
Searched
How are loads and forces calculated in the steady state analysis? (/home/mich/apps/CASEToolBox/casetoolbox/casestab/tests/DTU10MW)
基于对代码的分析，我可以确认是的，这个文件确实有载荷计算。让我详细说明：
载荷计算分析
1. 气动载荷计算
￼
# 从BEM数据中提取载荷信息bem_data = models[iops].rotors[0].create_data_for_BEM_results()# 面内力 (Inplane force) - 单位：kN/maxs1[1,1].plot(bem_data[:,0], bem_data[:,5]*1e-3, '.-',                hs2_ind[:,0], 1e-3*hs2_ind[:,6], '.-')# 面外力 (Out-of-plane force) - 单位：kN/m  axs1[2,1].plot(bem_data[:,0], bem_data[:,6]*1e-3, '.-',               hs2_ind[:,0], 1e-3*hs2_ind[:,7], '.-')
2. 载荷类型
A. 分布载荷
面内力 (Inplane force): 叶片平面内的气动载荷
面外力 (Out-of-plane force): 垂直于叶片平面的载荷
单位: kN/m (每米长度的力)
B. 整体载荷
￼
# 功率载荷pwr = rotor_models.pwr  # 功率 [W]# 推力载荷  thr = rotor_models.thr  # 推力 [N]
3. 载荷计算过程
A. 稳态计算
￼
# 执行稳态计算，包含载荷计算rotor_models.steady_state_computation()
B. 载荷数据提取
# 提取变形状态数据（包含载荷信息）
defl_data = models[iops].substructures[1].create_data_for_deflection_state()

# 提取BEM结果（包含气动载荷）
bem_data = models[iops].rotors[0].create_data_for_BEM_results()
5. 载荷输出
A. 图表输出
deflforce{:d}ms.png: 包含载荷分布图
显示面内力和面外力的沿叶片展向分布
B. 数据文件
# 保存稳态结果（包含载荷数据）
rotor_models.save_steady_state_results('/' + res_folder + '/ind/')




"""

pass

